invlogit <- function(xb) {
  1 / (1 + exp(-xb))
}